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Abstract 



Spike-timing-dependent plasticity (STDP) with asymmetric learning windows is 
commonly found in the brain and useful for a variety of spike-based computations such 
as input filtering and associative memory. A natural consequence of STDP is estab- 
lishment of causality in the sense that a neuron learns to fire with a lag after specific 
presynaptic neurons have fired. The effect of STDP on synchrony is elusive because 
spike synchrony implies unitary spike events of different neurons rather than a causal 
delayed relationship between neurons. We explore how synchrony can be facilitated by 
STDP in oscillator networks with a pacemaker. We show that STDP with asymmetric 
learning windows leads to self-organization of feedforward networks starting from the 
pacemaker. As a result, STDP drastically facihtates frequency synchrony. Even though 
differences in spike times are lessened as a result of synaptic plasticity, the finite time 
lag remains so that perfect spike synchrony is not realized. In contrast to traditional 
mechanisms of large-scale synchrony based on mutual interaction of coupled neurons, 
the route to synchrony discovered here is enslavement of downstream neurons by up- 
stream ones. Facihtation of such feedforward synchrony does not occur for STDP with 
symmetric learning windows. 

Keywords: spike-timing-dependent plasticity, synchronization, feedforward net- 
works, complex networks 
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1 Introduction 



In many neural circuits, synaptic plasticity depends on relative timing of presynaptic 
and postsynaptic spikes, which is known as spike-time dependent plasticity (STDP) 



(Gerstner et al., 1996 Bell et al., 1997 Markram et al., 1997 Bi and Poo, 1998 



Zhang et al., 1998). Specifically, long-term potentiation (LTP) ensues when a presy- 



naptic neuron fires slightly before a postsynaptic neuron (of the order of 10 ms), whereas 
long-term depression (LTD) is elicited in the opposite case (solid line in figured]). This 
STDP rule promotes causal relationship between a pair of neurons in the sense that the 
strength of a synapse that contributes to generation of postsynaptic spikes is reinforced. 
Computationally, STDP is useful for synaptic competition (Kempter et al., 1999 



Song et al., 2000[ |van Rossum et al., 2000[ |Song and Abbott, 2001]), coincidence de- 



tection ( Gerstner et al., 1996 ), spike-based associative memory (e.g. Lengyel et al. 
2005), implementation of the synfire chain ( Horn et al., 2000 Levy et al., 2001 ) 
generation of reproducible spatiotemporal spike patterns (Izhikevich et al., 2004 



|Izhikevich, 2006 ), selection of earlier inputs, to name a few (e.g. Gerstner and Kistler, 
2002). 

Related to neural computation, coincident firing of multiple neurons in the oscil- 
latory regime is found in many parts of the brain and believed to play an important 



role ( [Singer and Gray, 1995^ |Ritz and Sejnowski, 1997^ [Buzsaki and Draguhn, 2004[ ). 
One could imagine that real neural networks learn to synchronize spikes of different 
neurons by STDP-related synaptic plasticity, as suggested by some modeling studies. 
However, contribution of STDP to spike synchrony may be limited. For example, 
STDP can lead to division of a neural population into clusters in each of which neu- 
rons fire in spike synchrony ( Horn et al., 2000 Levy et al., 2001 ). This self-organizing 
process actually necessitates homogeneous synaptic transmission delays for different 
synapses and puts a strong restriction on the firing period. If there is just one cluster 



3 



(all neurons firing synchronously), the firing period has to be equal to the synaptic 
transmission delay. Similarly, if there are two clusters, the first cluster excites a syn- 
chronous volley in the second cluster after the synaptic transmission delay. The second 
cluster reexcites the first cluster in a similar way. The firing period is equal to the 
synaptic transmission delay multiplied by the number of clusters, which seems restric- 
tive. Alternatively, coincident firing is achieved via STDP if the amount of LTP and 
that of LTD that are caused by a presynaptic and postsynaptic spike pair are perfectly 
balanced ( Karbowski and Ermentrout, 2002 ). Evolution of coincident firing survives 
heterogeneity in neurons and in the amount of plasticity. However, how coincident 
firing is affected by the imbalance between LTP and LTD remains to be explored. 

Coincident firing in recurrent neural networks may not be established through 
STDP. In general, synchronous firing can be induced by sufficiently strong coupling be- 



tween elements ( [Kuramoto, 1984] |Pikovsky et al., 2001[ |Gerstner and Kistler, 2002| ). 
By contrast, STDP cannot strengthen the synaptic weights between two neurons bidi- 
rectionally. An increase of the synaptic weight in one direction implies a decrease in 
the opposite direction, and stability requires that the net decrease and the net increase 
are roughly balanced ( Song et al., 2000 Song and Abbott, 2001 ). Therefore, STDP 
does not necessarily enhance mutual interaction. Indeed, in recurrent neural networks, 
STDP does not necessarily support synchronous firing ( Masuda and Aihara, 2004 ). It 
rather reinforces reproducible spatiotemporal spike patterns composed of causal spike 
pairs of different neurons ( Izhikevich et al., 2004 Izhikevich, 2006 ). 

We examine possible mechanisms of STDP-induced synchrony in recurrent networks 
of oscillatory elements. We distinguish two types of synchrony using the terminology of 
coupled oscillators. One type is phase synchrony, which is equivalent to spike synchrony. 
When neurons are in phase synchrony, they share spike timing. The other weaker notion 
is frequency synchrony in which neurons possibly with different intrinsic firing rates 
share a common firing rate. Frequency synchrony does not imply phase synchrony. The 
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spike time of the postsynaptic neuron can differ from that of the presynaptic neuron. 
According to these definitions, the previous studies cited above, which relate STDP to 
synchrony, regard phase synchrony. 

STDP may be more relevant to frequency synchrony. For example, STDP pro- 
motes frequency synchrony, but not phase synchrony, in a hybrid circuit of an aplysia 



abdominal neuron and an emulated neuron ( Nowotny et al., 2003 ). Unidirectional con- 
nectivity from the emulated neuron to the aplysia neuron eventually forms. Numeri- 



cal simulations of two coupled Ho dgkin- Huxley neurons (Zhigulin et al., 2003) and of 



large neural networks ( [Zhigulin and Rabinovich, 2004] ) also support the notion that 
frequency synchrony is facilitated by STDP. 

In the present work, we show that the standard STDP facilitates frequency syn- 
chrony to a great extent, particularly when LTP and LTD are roughly balanced. To 
examine how heterogeneous neurons interact to produce possible synchronization, we 
analyze networks of oscillators with a pacemaker. The pacemaker has a distinct nat- 
ural frequency, is not affected by other oscillators, and sets the rhythm to influence 
other neurons. No matter whether a pacemaker is realized by a network or a sin- 
gle neuron, existence of pacemaker neurons is suggested in, for example, the basal 



ganglia (Plenz and Kitai, 1999) and respiratory networks in the pre-Botzinger com- 



plex (Ramirez et al., 2004). Furthermore, many neurons (Hutcheon and Yarom, 2000) 



and recurrent microcircuits ( Jefferys et al., 1996 ) are intrinsically oscillatory (also see 
e.g. Singer and Gray, 1995). and their rhythmic activities are resistant to perturbation. 
These neural networks and single neurons can also serve as pacemakers. With frozen 
and sufficiently strong synapses, oscillator networks with pacemakers allow frequency 



synchrony ( Kori and Mikhailov, 2004 ). We show that STDP considerably facilitates 
frequency synchrony of pacemaker systems by establishing feedforward network struc- 
ture whose root is the pacemaker. 

For analytical tractability, we mostly deal with networks of phase oscillators in 
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which couphng strength evolves according to STDP. Coupled phase oscillators approx- 
imate various natural systems composed of self-sustained oscillators with weak coupling 



(Winfree, 1980 Kuramoto, 1984 Glass and Mackey, 1988), including pulse-coupled 



neurons (Kuramoto, 1991 Hansel et al., 1993 Hansel et al., 1995 Kori, 2003). We 



introduce the model in section [2] and analyze simple cases of two connected oscillators 
in section [3l We numerically analyze larger oscillator networks with STDP in section HI 
In section [5l we numerically simulate pulse-coupled pyramidal neuron models to show 
that our results obtained for coupled phase oscillators qualitatively apply to spiking 
neuron models. 

2 Model 

We analyze a network of n phase oscillators. One oscillator is assumed not to be dis- 
turbed by the other n — 1 oscillators. We designate this special oscillator as pacemaker 
and use the term oscillator to refer to the other n — 1 elements. The pacemaker has 
natural frequency Vt and phase 0o ^ [0, 27r). The other oscillators are assumed to have 
the identical natural frequency w, and the phase of the i-th oscillator {1 < i < n — 1) 
is denoted by 0j G [0, 27r). We identify 0i = and (pi = 2-n (0 < i < n — 1). We write 
E E ii there is a synaptic connection from oscillator i to oscillator j. In other 
words, E is the set of edges of the underlying neural network. As in real neural net- 
works, connectivity is asymmetric in general so that [i-,]) G E does not imply (j, i) G E. 
A pair of connected oscillators interact via sinusoidal coupling, which usually promotes 



synchrony (Kuramoto, 1984) 



Dynamics for fixed synaptic strengths are represented by: 



a (11 



k = ^^ + 7^: fi'jiSin(0j - {l<i<n-l) (2) 

where (k) is the average number of incoming edges per oscillator. The coupling strength 
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Qji is associated with synapse (j, i). We note that Qio, which is the synaptic weight from 
an oscillator to the pacemaker, does not affect the network dynamics: the pacemaker 
is not perturbed by external input. However, we will monitor QiQ to examine how this 
connection evolves as synaptic plasticity goes on. 

We assume that Q > u unless otherwise stated. By rescaling the timescale and the 
coupling strength, we set Q = u+l without losing generality. To set the values of Q and 
u, we take care of two subtle factors. First, a small u would yield backward rotation 
by the effect of coupling. This is because the second term of the right-hand side of 
equation ([2]) can be large negative to overwhelm the first term. Then, the condition 
(pi < may be satisfied for long enough time to elicit backward firing. This is unrealistic 
as a neuron model. Second, we avoid a pair of Q and u that accommodates the relation 
MiQ = M2UJ with small integers Mi and M2(Mi ^ M2). In such a situation, resonant 
behavior appears when the pacemaker and the oscillators are decoupled through STDP 
and has a pathological effect (see the explanation after equation flTB]) for more details). 
The resonant firing is ruled out by dynamical noise in many real neural networks. 
However, we have to carefully specify Q and u in the present work because we do not 
assume noise for analytical tractability. Keeping these caveats in mind, we set i7 = 9.1 
and u = 8.1. 

Spike time is defined to be the time when the 0j crosses 0. Synaptic update based 
on STDP takes place based on a pair of nearest presynaptic and postsynaptic spike 
times, without paying attention to remote spike pairs (see arguments in e.g. Froemke 
and Dan, 2002). We compare the upshot of two types of STDP rules for synapse gji 
€ E), namely, asymmetric STDP and symmetric STDP. 

Asymmetric STDP is modeled as follows. LTP is induced if a presynaptic firing 
(spike of oscillator j) precedes a postsynaptic firing (spike of oscillator i). In the 
opposite case, LTD occurs. We denote the presynaptic (postsynaptic) spike time by 
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tpre {tpost)- A spike-pair event modifies the synaptic weight: gji — )■ gji + Agji, where 
Ag,, = ( ^''P (7^ 7f^2 \ ^'"-^ ^ (3) 

I A. exp ( ^T"^ j , tpj-e ^ tposti 

under the hmitation gji G [0, gmax]- A sample learning window is indicated by the solid 
line in figure [H The width of the learning window is specified by r, which is known 



to be of the order of 10-20 ms (Bi and Poo, 1998 Zhang et al., 1998). We confine 



ourselves to the regime in which firing rates are not very large (5-20 Hz), as is true 
for many pyramidal neurons. Then, r is several times smaller than the characteristic 
interspike interval T = 50-200 ms. We thus set 

1 27r T 

r = - X — = -. 4 

6 n 6 ^ ^ 

For completeness, we assume that tpre = tpost does not induce plasticity. 

We assume that synaptic weights evolve so slowly that we can solve equation ([2]) by 
regarding the synaptic weights as constant. This assumption is valid if A'^Q, A^Q ^ 
g'^ for the following reason. Because the relative phase relationship determines the 
evolution of synaptic weights, we should compare the typical timescale of the relative 
phase dynamics with that of synaptic plasticity. The former is the inverse of typical 
synaptic weight gQ. By introducing dimensionless synaptic weight g/go, we find from 
equation that the timescale of dimensionless synaptic plasticity is the inverses of 
A^Q/gQ and A^Q/gQ. The two timescales are separated if A'^Q/g^, A"Q/gQ ^ g^, 
which leads to A^Q, A^Q <^ g'^ = g"^. On the slow timescale of synaptic plasticity, we 
can set A^ = 1 by rescaling the time, so that only the ratio A'^ /A^ is relevant. For the 
stability, A'^ /A^ must be balanced. This ratio is assumed to be slightly smaller than 
unity according to previous literature dSong et al.,^000l [Song and Abbott, 2001[). 



Most of our theoretical efforts are invested in asymmetric STDP because many 
pyramidal neurons show asymmetric STDP. However, symmetric STDP, in which 
the synaptic update rule depends only on \tpre — tpost\, is also found in some experi- 
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ments. Particularly, the learning window is often shaped like a mexican hat in excita- 



tory synapses; small (large) \tpre — tpost\ induces LTP (LTD) ( |Nishiyama et al., 2000 



Abbott and Nelson, 2000} [Shouval et al., 2002]) . Symmetric learning windows have 



been found for inhibitory synapses (Woodin et al., 2003) and for the amount of LTD 



in excitatory synapses (Dan and Poo, 1992). We numerically analyze networks with 



symmetric STDP in section HI We adopt superposition of two gaussian distributions as 
the symmetric learning window, as depicted by the dotted line in figured] A spike-pair 
event modifies the synaptic weight: gji — )■ gji + Agji, where 

\ f (tpre tpost) \ ^ / {tpre tpost) \ /-\ 

= 75^ I ^^^j - [ ^^^j ■ 

with cr+ < a^. We set = 0.6r and a~ = 2a^ = 1.2t so that the timescale of the 
symmetric learning window is comparable to that of the asymmetric learning window 
defined in equation The values of A'^ and A" are assumed to be the same as those 
for asymmetric STDP so that gji is bounded. 

3 Analysis of Small Networks with Asymmetric 
STDP 

We begin with small networks of two oscillators with asymmetric STDP. For these 
networks, how much initial coupling is necessary for synchrony can be analytically 
evaluated. 

3.1 One Pacemaker and One Oscillator 

We deal with the case n = 2, namely, a network of one pacemaker and one oscillator. 
Because the connection from the oscillator to the pacemaker does not affect the dy- 
namics of the pacemaker, it suffices to consider the unidirectional case. The network 
is schematically shown in figure HJ^a). We write g = g^i to simplify the notation. The 



9 



short-term dynamics in which g is regarded to be constant are described by 



00 = ^, (6) 

01 = uj + gsm{(j)o- (pi). (7) 

With ip = (pQ — (pi, equations and ([7]) reduce to 

ip = Q — u — gsinip. (8) 

Based on the assumption that synaptic plasticity occurs much more slowly than firing, 
we perform quasistatic analysis. For a frozen g, let us derive the average angular 
frequency of the oscillator denoted by a). If g > Q — u, the pacemaker and the oscillator 
are in frequency synchrony, i.e. u = Q. lfO<g<Q — co, equation (E]) is equivalent to 

dip 



I 



fl — uj — gsinip 
Integration of equation ([9]) over a cycle yields 



dt. (9) 



2tt fT , f^^ dip 

' dt ' 



Q — oj Jo Jo Q — oj — g sin ip 

27r 



(10) 



which results in 



Co = n- ^{n-ooy - g^. (11) 

Note that u < < Q. 

The direction and the amount of synaptic plasticity induced by a single spike-pair 
event is determined by tpost — tpre- We estimate tpost — tpre in terms of ip as follows. 
Suppose that the phase difference is equal to ip when the pacemaker fires. Then, it 
approximately takes tpost — tpre = 4'/^ fo^' the oscillator to fire. In this case, LTP 
is induced because the pacemaker is presjTiaptic to the oscillator. The pacemaker 
and the oscillator can fire in the opposite order. If the phase difference is ip when 
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the oscillator fires prior to the pacemaker does, the pacemaker spends approximately 
tpost — tpre = before firing. In this case, LTD is induced. Because we confine 
ourselves to the case in which Q does not deviate so much from u, we approximate 
tpost — tpre = ip/fl regardless of the order of firing. 

In equation (jl]), we assumed that the decay rate of the learning window r is suffi- 
ciently smaller than T/2, which corresponds to phase vr. Therefore, the amount of LTP 
is negligible for tpost — tpre — tt/o; or larger, and the LTP rule is effective only when 
< ip < TT. By the same token, the LTD rule is effective only for — vr < ip < 0. Using 
these approximations, we aim to describe the dynamics of synaptic plasticity in terms 
of phase variables. The amount of plasticity given by equation ([3]) can be rewritten as 

, A+ exp f-^) , < ^ < vr, 01 = 
Ag= { V < 12 

^ -A- exp (^) , -TT < ^ < 0, 00 = 



where 0i = and 0o = indicate the postsynaptic spike time for an LTP event (the 
spike time of the oscillator) and that of an LTD event (the spike time of the pacemaker), 
respectively. 

We denote by g{0) the initial synaptic weight. If g{0) > Q — u, fast dynamics have 
two steady states given by 

ip* = arcsin (^"^^j • (13) 

The solution with 7r/2 < -i/;* < vr is unstable, and hence the fast dynamics converges to 
ip* satisfying < ip* < tt /2. Therefore, STDP induces potentiation of g. Then, ip* for 
an altered g becomes even smaller, which induces further potentiation of g. Eventually, 
9 = Qmax is achieved. In sum, if 

giO)>gc = n-cu, (14) 

the pacemaker and the oscillator will synchronize quickly without plasticity. The STDP 
does not break frequency synchrony. Note that STDP generally decreases the phase 
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difference ip, but does not tend to (no pliase synclirony) unless Q = u. Alternatively, 
if gmax = CO, g diverges, and ip goes to 0. 

Does asymmetric STDP facilitate frequency synchrony? If g{0) < g^, the oscillator 
is initially not entrained by the pacemaker. Then, ip slips. By averaging over many 
spike-pair events, we represent the synaptic dynamics on a slow timescale by 

Jt such that o<i/'<7r Jt such that -7r<i/i<o 



/TT 



;i5) 



0. re'^/"-f - (16) 

Jo \gc-giiini) gc + gsmipj 

Derivation of equation fllSp requires the nonresonant situation. If MiQ = M2UJ holds 

with small Mi and M2, the dynamics become periodic with a rather small period when 

the oscillator decouples from the pacemaker due to STDP. If this were the case, the 

dependence on the initial condition does not vanish permanently. In other words, ip 

conditioned by 0o = or 0i = in equation (fT2i) would take only limited values. Then 

the distribution of conditioned by a spike event would deviate from the unconditioned 

distribution of ip. With our choice ofQ = 9.1 and u = 8.1, the effect of such a resonance 

is very small. 

In the region of g where g > holds, the RHS of equation (fT6i) increases monotoni- 
cally with g. If g{0) is greater than the value that makes the RHS equal to zero, which 
we denote by gc-stdp, we obtain g > 0. Under this condition, g continues to increase, 
and ip* decreases. This makes g in equation f|T6|l even larger. This positive feedback 
lasts until g > gc is eventually satisfied. As a result, frequency synchrony is elicited 
by STDP. However, if g{0) < gc-stdp, 9 converges to the lower bound 0, so that the 
oscillator is disconnected from the pacemaker. 

We bound gc-stdp as follows: 



RHS of equation f lTB]) 



g sin 'il) 



^ gc ^ \ 9c 
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9c Jo y ^ + fj ^0 

{A- - A+)Qt{1 - e-"/^0 g A+ + A- Q'^t\1 + e-"/^") 



9c 9l i + i 1 + 



(17) 



The value of g that makes the RHS of the above equation zero gives an upper bound 
of gc-stdp- When y4+ and A~ are balanced (^4+ = A~), we obtain 

^ Q.(l+e-V»^) ^ A--A+(l + l^V)(l-e-/"-) 

^c-stdp _ _ ^-_^+ (i+n2r2)(i-e-V"^) ~ A- + firfl + e-^/f^^) ^ ^ 

Because fir is assumed to be of the order of n (see equation (|1])), Qr = 0(1). In 
addition, when ^4+ = A~, the inequalities in equations f ll7p and f llSp nearly hold with 
the equalities. In such 

gc-stdp oc ^1 - (19) 

which implies that gc-stdp is much smaller than gc when A'^ = A~ . Particularly, gc-stdp 
is extinguished when A'^ > A~ . 

In figure [3l we plot gc-stdp evaluated by numerical integration of equation f|T6|) and 
the approximation given by equation (fTSi) . We also plot gc-stdp obtained by numerical 
simulations of our model (equations ([3]), ([S]), and ([7])), in which gc-stdp is determined by 
varying the initial synaptic weight g{0). The evaluation by equation flTBl) (solid line) is 
in good agreement with gc-stdp obtained by numerical simulations of the model (circles) 
for a broad range of A'^ /A~ . As expected, the approximate estimation by equation flTS]) 
(dotted line) also agrees with the numerical data (circles) when A'^ /A^ is close to unity. 
In conclusion, asymmetric STDP drastically enhances frequency synchrony. 

Regarding symmetric STDP, for values of g such that ip falls in the positive learning 
window (refer to the dotted line in figured]), g is strengthened to eventually exceed gc- 
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Therefore, frequency synchrony is facihtated. To what extent synchrony is promoted 
depends on the width of the learning window. 

When Q < u, there are two solutions ip* G (— 7r,0), one of which is stable. Then, 
STDP elicits LTD. Even though ip* changes, the relation — tt < ^* < is preserved until 
the pacemaker and the oscillator get disconnected. As a result, frequency synchrony 
does not happen. 

3.2 Two Oscillators 

To examine how the connectivity between a pair of oscillators evolves in a large network, 
we analyze the following toy model of two bidirectionally coupled oscillators: 

4>i = u + Au + gism{(j)2 - (f)i), 

h = a; + 5f2sin(0i - 02), (20) 

where gi, g2 G [0,gmax]- The network is depicted in figure |5|^b). Now two oscillators 
influence each other, which contrasts to the case of the pacemaker-oscillator network 
examined in section 13.11 The term Ao; represents the mismatch in natural frequen- 
cies. Although the oscillators are identical in our original model (see equation ([2])), 
we introduce Au because of the following reason. In oscillator networks with a pace- 
maker, the oscillators are not completely phase synchronized. The oscillators directly 
connected to the pacemaker are the first to fire after the pacemaker does. Then, other 
oscillators adjacent to those connected to the pacemaker fire after some delay, and so 
forth. Therefore, the oscillators closer to the pacemaker tend to have more advanced 
phases, and the distribution of the phases is associated with the hierarchical organiza- 
tion of the network. Imagine two oscillators coupled unidirectionally or bidirectionally 
in a large network. We denote one that fires first and the other by oscillators 1 and 2 
respectively. Precisely, the difference in the firing timing stems from complex effects of 
coupling with other oscillators. For analytical tractability, here we replace such effects 
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by the frequency mismatch Au, by which the difference in the firing timing can be 
easily introduced. As shown in the following, for Au > 0, oscillator 1 tends to fire 
in advance of oscillator 2. Accordingly, we regard that oscillator 1 is closer to the 
pacemaker than oscillator 2. 

We analyze the model given by equation fl20l) . By introducing ■?/' = 0i — (/)2, we 
derive 

ij = Aoj - {gi + g2)sm'ip. (21) 

If 

gi{0)+g2iO)>Au, (22) 
two oscillators are locked with phase lag < ip* < 7r/2, where 

tJj* = arcsin f — — | . (23) 

\9i + g2j 

If synaptic plasticity is absent, equation fl22]) gives the condition for frequency syn- 
chrony. 

In contrast to the network of one pacemaker and one oscillator analyzed in sec- 
tion I3.lt equation ( 12^ does not guarantee that frequency synchrony is maintained 
throughout STDP. When equation fl22|) is satisfied, the synaptic dynamics are written 
as 

g, = -A-exp(^-^y (24) 
g2 = A+exp(^-^y (25) 

where 

UJ = UJ -\ 26 

91+92 

is the frequency common to the two oscillators. Because A'^ < A~ , gi + g2 decreases 
with time. The oscillators desynchronize in frequency if gi + g2 > Au is violated 
via synaptic plasticity. For sufficiently small gi{0) + (72(0), the two oscillators are 
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disconnected even from the beginning. In these cases, ip shps due to the absence of 
frequency synchrony. 

The average frequencies of the two oscillators out of frequency synchrony are cal- 
culated as 



UJ2 



+ 92)' 



91 + 92 



g2Auj - g2jAuj^ - {91 + 92Y 



The synaptic weights evolve according to 



9i 



1 

2^ JO 
1 



91 + 92 



(27) 
(28) 



Auj + {gi + g2) sin ^ Au - {gi + g2) sin ^ 

A+ A- 



92 



27r Jo 

2tt Jo 
1 



Au + {gi + 5(2) sin ip Auj - {gi + 5(2) sin ^p 

' dip 



dip, (29) 



271 Jo 



Auj - {gi + g2) sin ip Auj + {gi + g2) sinip 



dip, (30) 



Auj - {gi+ g2) sin ip Auj + {gi + g2) sin ip _ 
where we approximated e"'^/'^^^ i^y g-iA/i^ir^ ^-^y section IXTl 

Since < is always satisfied, gi eventually reaches 0; backward connectivity from 
a downstream oscillator to an upstream oscillator is eliminated. Whether a 'forward' 
connectivity from the upstream oscillator to the downstream oscillator survives relies 
on g2{i) where i is the time gi reaches 0. If g2{i) is larger than gcstdp obtained in 
section 13.11 frequency synchrony will be eventually established. In this case, the final 
oscillation frequency is equal to u + Au so that oscillator 2 is enslaved by oscillator 1. 
If 5^2 (^) < gcstdp, the two oscillators are finally disconnected. 

The critical value of (72(0) above which frequency synchrony occurs is plotted in 
figurelUfor different values of A'^ /A~ and gi{0). The critical (72(0) decreases with gi{0), 
implying that a large gi (0) enhances frequency synchrony. Such backward connectivity 
transiently serves to keep ip small so that forward connectivity is enhanced. However, 
only the synapse from the faster oscillator to the slower oscillator survives eventually. 
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The feedforward network is not created via symmetric STDP, by which gi and 
g2 evolve in the same direction. When initial mutual connectivity is strong enough, 
synchrony is established so that the two synaptic weights are saturated {gi = g2 = 
Qmax)- Then, based on equation (l26ll . the common firing frequency is equal to (X'+Au;/2, 
but not to the frequency of the faster oscillator (= a; + Aw). 

The effect of LTP-LTD balance is also shown in figure |H When is close to A~ , 
critical (72(0) is lowered, and entrainment occurs easily. Even when LTP is rather weak 
compared to LTD {A'^ /A~ = 0.8, thinnest line), the critical g2 is much reduced from 

gc- 

4 Numerical Results for Large Networks 

So far, we have analyzed small networks composed of two elements only. In this section, 
we examine how frequency synchrony can be facilitated by STDP in larger networks. 
In particular, we compare gc and gc-stdp and also investigate evolution of network 
structure. To this end, we numerically simulate randomly connected n = 100 elements 
(99 oscillators and one pacemaker) based on equations ([3]), ([1]), and ([7]). 

4.1 Initial Setup 

We generate a directed random network as follows. Starting from a set of n isolated 
vertices, we add a directed edge that connects a randomly chosen pair of oscillators. 
We forbid multiple directed edges between a pair of oscillators and self loops, i.e. edges 
whose origin and destination are identical. This procedure is repeated n (k) times. In 
the following, we assume that (A;) = 10. In other words, each oscillator is presynaptic 
to 10 other oscillators and postsynaptic to the same number of oscillators on average. 

We define the distance /j of oscillator i from the pacemaker by the smallest number 
of directed edges necessary to reach from the pacemaker to oscillator i. For example, 
the number of the oscillators at distance 1 is equal to those that receive direct synaptic 
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contacts from the pacemaker. Therefore, about (k) oscillators have distance 1. Among 
the other oscillators, those receiving an edge from an oscillator with distance 1 have 
distance 2. The depth L of a network is defined as the distance averaged over all the 
oscillators: L = Y17=i U/in — !)• 

The initial phases 0i (0 < ? < n — 1) are taken independently from the uniform 
density on [0,27r). For all the synapses, we initially set gji = g{0). 

For a specific random network used in the following numerical simulations, we 
obtained L = 3.22. For this network, we numerically found that frequency synchrony 
happens without synaptic plasticity if g > gc — 100.7. In this case, the pacemaker first 
fires in each cycle, and oscillators with smaller distances tend to fire with smaller lag 



with respect to the pacemaker (Kori and Mikhailov, 2004) 



4.2 Measured Quantities 

We define the degree of frequency synchrony r = r{[ti,t2]) for a time interval [^1,^2]- 
The mean frequency of each oscillator for this time interval is equal to 



We note that ujq = Q. Then, the synchrony measure is defined by 



(32) 



— tu 

When the oscillators are in frequency synchrony with the pacemaker, the mean fre- 
quency of the oscillators J2i=i ~ 1) is equal to the frequency of the pacemaker Q, 
and we have r = 1. If the pacemaker does not at all affect the other n — 1 oscillators, 
the oscillators fire at their natural frequency u, and we have r = 0. We divide the 
total simulation time into consecutive bins of the width t2 — ti = 100 for oscillator 



simulations in sections 14.31 and 14.41 

More microscopically, we inspect possible formation of feedforward chains originat- 
ing from the pacemaker. To quantify this, we track several order parameters derived 
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from synaptic weights. The first is the depth L extended to networks with heteroge- 
neous synaptic weights in the following way. Let us consider a path from the pace- 
maker to oscillator i. A path is equivalent to a chain of directed synapses: {jo,ji) G E, 
(ji, js) e E, ijk,-i,jki) e E, where jo = and jk, = i- The length of this path 
is given by Y!k=o Qmax/ gj^jk+i- The distance U of oscillator i from the pacemaker is 
the shortest path length among all possible paths from the pacemaker to oscillator 



i ( Braunstein et al., 2003 ). This definition generalizes the prior definition for networks 
with unit synaptic weights. The redefined distance is associated with how much a 
downstream oscillator is influenced by the pacemaker. The depth of the network is 
again defined by L = Y^^Zi k/{n — 1) and measures effective proximity of the oscilla- 
tors from the pacemaker. By definition, the generalized L is equal to or larger than 
L of the unweighted network, with equality realized only when gji = Qmax for all the 
synapses that appear in the shortest paths. 

A synaptic connection (j, i) with Ij < li [Ij > li) is a forward (backward) connection 
in the meaning that it complies with the feedforward chain emanating from the pace- 
maker. Accordingly, we define the amount of forward connection Wf, that of backward 
connection Wf,, and that of lateral connection wi by 

_9ji 

1 

-e<.li — lj <e 



Gi = E (35) 



Summation is taken over the pairs of oscillators forming synapses ((j, i) ^ E). Note that 
Gi quantifies the connection between oscillators whose distances from the pacemaker 
are approximately equal. The number of synapses in the network (= n{k)) gives 
normalization, and thus < G/, G^, Gp < Qmax- The average synaptic weight is given 
by < Gf + Gh + Gp < Qmax- The tolerance level is chosen to be sufficiently small: 
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e = 0.05. We also define local quantities to evaluate feedforwardness. The average 
weight of the synapses postsynaptic to the pacemaker is denoted by Gj. This is equal 
to the average of goi over i with {0,i) G E. This corresponds to g used in section [STTl 
Similarly, the average weight of the synapses presynaptic to the pacemaker is denoted 
by G^. This is equal to the average of gio over i with (z, 0) G E. We note that 

^0 



0<G<},Gt<gr, 



4.3 Asymmetric Learning Window 



We apply asymmetric STDP with LTD slightly stronger than LTP (Song et al., 2000 
Song and Abbott, 2001^ : A+ = 0.009 and = 0.01. By setting g^ax = 15 < gc, ho- 
mogeneous enhancement of all the synapses does not lead to synchrony. We determine 
Qc-stdp by running numerical simulations with various values of the initial synaptic 
weight gji = g{Q). 

Because A'^ is close to A~ , our results in section [3] predict the following. 



• As shown in section 13.21 backward connection will be eliminated via the asym- 
metric STDP so that and G° decrease. 

• The unidirectional connection between the pacemaker and the oscillator will be 
easily established (section [XT]) . As a result, a feedforward chain rooted at the 
pacemaker is expected to form. 

• gc—stdp IS much smaller than gc- 

Dynamics of synaptic-weight order parameters for g{0) = 0.7 are shown in fig- 
ure [5]^ a, b). The average synaptic weight (dotted line in figure [5](a)) increases in the 
initial stage because some synapses between the oscillators are potentiated. However, 
its stationary value is much smaller than the maximal possible value {gmax = 15). There 
is no selective potentiation of forward synapses (G/, thick solid line in figure [5]^a)) or 
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depression of backward synapses {Gh, thin solid line in figure El^a)). The forward con- 
nectivity from the pacemaker (Gj, thick line in figure El^b)) also degrades with time. 
Eventually, the oscillators disconnect from the pacemaker, which is observed as indef- 
initely growing L (uppermost line in figure EJ^e)). Accordingly, frequency synchrony 
between the pacemaker and the oscillators is not achieved; figure [5]^f ) indicates that r 
stays near 0. 

By contrast, frequency synchrony without phase synchrony is established when 
g{0) = 1.5, as supported by the rastergrams in figures [6]( a) and[6]^b) corresponding to 
initial and final periods of simulations, respectively. More in detail, forward connec- 
tivity Gf (thick solid line in figure |5]^c)) and (thick line in figure [5]^d)) grow toward 
dmax to result in frequency synchrony at t = 12500 (figure EJ^f)), accompanied by a 
decrease in L (lowermost line in figure |5]^e)). Backward synapses directly projecting 
to the pacemaker are pruned in an initial stage (GjJ; thin line in figure [5]^d)). It takes 
longer time for Gh to decay (thin solid line in figure EJ^c)). Although randomness in 
the initial condition blurs the phase transition, we estimate Qc-stdp — 0.9 based on 
figures [S](e) and (f). Consistent with the results in section [XT| Qc-stdp is much smaller 
than = 100.7. 

Frequency synchrony is made possible by combined effects of sufficiently large for- 
ward weights and sufficiently small backward weights. It is not an immediate conse- 
quence of increased average synaptic weights; achieving synchrony merely by homo- 
geneously strong synapses necessitates G f + Gi, + Gp > Qc- Because of our choice of 
Qmax (< Qcli homogeneous LTP does not induce frequency synchrony. Elimination of 
backward weights is essential for frequency synchrony. The final network structure re- 
constructed from synapses with gji > g{0) = 1.5 is shown in figure |5]^g) , with forward 
and backward edges shown by thin and thick lines, respectively. Few backward edges 
survive asymmetric STDP. The network is close to a feedforward network rooted at 
the pacemaker, which enslaves the oscillators. 
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We remark that detailed behavior of the network order parameters varies according 
to the value of e and the definition of the distance /j, which is inherently arbitrary for 
weighted networks. However, the general tendency that forward synapses are poten- 
tiated and backward synapses are depressed for g > gc~stdp is observed irrespective of 
these details. For a reasonably defined U, whether L increases or decreases (figure [5](e)) 
and whether feedforward networks form (figure [5](g)) are determined independently of 
the definition of k. 

We have performed additional numerical simulations in which every presynaptic 
spike spends for r/5 = T/30 before exciting the postsynaptic neurons. Because T = 
50-200 ms, the corresponding synaptic delay is equal to a few milliseconds. The value 
oi Qc-stdp hardly changes with this synaptic delay (results not shown). 

To examine the effect of heterogeneity in oscillators, we pick the intrinsic frequency 
of each oscillator from the gaussian distribution with mean w = 8.1 and standard devi- 
ation 0.1. Time courses of r are shown in figure El^h) . We estimate 1.2 < Qc-stdp < 1-5, 
implying the robustness of our results against heterogeneity. We remark that r con- 
verges to a positive level when g < Qc-stdp- This is because, even if the oscillators 
disconnect from the pacemaker, some oscillators form feedforward networks of small 
size in which fast oscillators entrain and speed up slow oscillators. If the heterogeneity 
is even larger so that some oscillators are as fast as the pacemaker, frequency syn- 
chrony seeding from the pacemaker would be difficult because the pacemaker and fast 
oscillators compete in entraining slow oscillators. 

4.4 Symmetric Learning Window 

Now we examine symmetric STDP. We set gmax = 200 > gc = 100.7 so that frequency 
synchrony with small phase lags is achieved if gji = gmax for all (j, i) G E. For the 
network same as that used in section 14. 3[ evolution of synaptic weights are summarized 
in figures [D^a, b) and (c, d) for g{0) = 140 and g{0) = 150, respectively. Because 
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the oscillators share an identical natural frequency, the phases are fairly close among 
them even with weak coupling. The synapses among these oscillators are potentiated. 
Accordingly, for both values of g{0), the average synaptic weight initially increases 
(dotted lines in figure [T]^ a, c)). Note that the forward weights {Gf, thick solid lines 
in figure [Tl^a, c)) and the backward weights {Gb, thin solid lines in figure [Tl^a, c)) 
are equally potentiated. Accordingly, the distance between the oscillators is initially 
shortened to result in a decrease in L (figure [7](e) , t < 2000). 

When g{0) = 140, the average synaptic weight stops increasing at a value slightly 
smaller than Qmax = 200 (dotted line in figure [7|( a)). This is because the synapses link- 
ing the pacemaker to the oscillators have not been potentiated. Actually, (thick line 
in figure [TJ^b)) and G^ (thin line in figure [7](b)) decrease to eventually decouple the os- 
cillators from the pacemaker. Consequently, L diverges (uppermost line in figure [Tl^e)), 
and frequency synchrony is eventually lost (figure [Tl^f)). 

When g{0) = 150, (thick line in figure [7|^d)) and G° (thin line in figure [Tl^d)) as 
well as Gf (thick solid line in figure [7](c)) and Gb (thin solid line in figure [7](c)) increase. 
Consequently, L continues to decrease to reach the minimum possible value for which 
gji = Qmax is achieved for most synapses (lowermost line in figure [7](e)). The synchrony 
measure r stays near unity throughout (figure [7](f)). In fact, approximate phase syn- 
chrony as well as frequency synchrony has been realized quickly. The synchrony arises 
not owing to STDP but to sufficiently strong initial coupling. 

Based on figures [Tl^e) and (f), which show time courses of L and r for several values 
of g{0), we estimate 145 < Qc-stdp < 146. This value of Qc-stdp is much larger than 
the case of asymmetric STDP and comparable to Qc = 100.7, that is, the threshold for 
frozen synapses. 

The fact that symmetric STDP does not promote frequency synchrony man- 
ifests the importance of the feedforwardness of networks. In general, forward 
synapses promote frequency synchrony, whereas backward synapses hamper it 
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( Kori and Mikhailov, 2004 ). Symmetric STDP does not independently control for- 
ward synaptic weights and backward synaptic weights. Consequently, it cannot get rid 
of backward synapses without sacrificing forward synapses. Final network structure 
is shown for g{0) = 150 in figure [Tl^g) . In contrast to the case of asymmetric STDP 
(figure [St^g) ) , many backward synapses (thick lines) remain. Under symmetric STDP, 
frequency synchrony is due to strong mutual interaction but not to formation of a 
feedforward network. 

5 Networks of Pulse-coupled Spiking Neurons 

To inspect whether the results obtained for coupled phase oscillators qualitatively ap- 
ply to more general neuron models, we numerically simulate pulse-coupled spiking 



neurons under STDP. We adopt a two-dimensional neuron model (Izhikevich, 2003 



Izhikevich et al., 2004). The subthreshold dynamics of the i-th. neuron are described 



by 

Vi = OMvf + 5Vi + 140 - Ui - Isyn,i - Iext,i, (36) 

iii = a{bvi-Ui), (37) 

where Vi is the membrane potential (mV), Ui denotes the recovery variable that evolves 
slowly relative to Vi, and the time unit is millisecond. The spiking mechanism is im- 
plemented by resetting the dynamical variables to (f m,) = (c, d) as soon as Vi exceeds 
30 mV. We set a = 0.02, b = 0.2, c = —65, and d = 8, which are standard parameter 



values for modeling pyramidal neurons (Izhikevich, 2003 Izhikevich et al., 2004). 



The input Iext,i and Isyn,i are the external bias input and the synaptic input, re- 
spectively. We set Iext,o = 8.4 and Iext,i = 8{l<i<n — 1). The inherent firing rate 
of the pacemaker (= 18.8 Hz) is about 5 % higher than that of the oscillators (= 17.9 
Hz). In figure [8], example traces of the membrane potentials of the pacemaker (solid 
line) and that of an oscillator (dashed line) are shown. 
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The synaptic input Isyn,i is composed of superposition of incident spikes from the 
neurons presynaptic to the i-th neuron. A presynaptic spike of the j-th neuron ((j, i) G 
E) is assumed to change the postsynaptic membrane potential according to the time 
course given by the alpha function: 

gj,s{t) = g,iaHe-''\ t > 0, (38) 

where t = corresponds to the spike time. We set a = 1, so that the unit synaptic 

input s(t) peaks at t = 1/a = 1 ms and then decays slowly. We set = 0.09, 

A" = 0.1, and r = 10 ms. 

The random network with n = 100 used in the following simulations are the same 

as that used in section HI We set the initial synaptic weight gji = g{0) for all the 

synapses. The initial values of Vi and Ui are independently chosen according to the 

uniform distributions on [—75, —50] and [—8, —6], respectively. Under these conditions, 

we track time courses of Gf, Gb, Gi, Gj, G^, L defined in section and r = r([ti, ^2]) 

redefined based on spike counts: 

J27=i (number of spikes from the i-th neuron) 
(number of spikes from the pacemaker) ^ ^ 

Frequency synchrony yields r = 1. If the pacemaker and the oscillators fire indepen- 
dently, r is the ratio of the single-neuron firing rate to the pacemaker firing rate. We 
set t2-ti = 10000 ms. 

With these parameter values, we first determined gc without STDP. We numerically 
obtained g^ = 1.6 for the network of one pacemaker and one oscillator, and gc = 38 
for the random network. In the following simulations with asymmetric STDP, we set 
9max = 35 < so that uniform increases in gji do not cause synchronization. Frequency 
synchrony requires feedforward network structure. 

For homogeneous initial synaptic weights g{0) = 5 and g{0) = 10, evolution of 
synaptic weights via asymmetric STDP is shown in figures [9]^a, b) and M^c, d), re- 
spectively. For g{0) = 5, Gf (thick solid line in figure M^a)) and G'j (thick line in 
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figure in^b)) do not grow during the course of plasticity, similar to figures [S](a, b). The 
oscillators disconnect from the pacemaker, and frequency synchrony is not realized 
(figure int^f)). For g{0) = 10, the forward connection from the pacemaker to the set of 
oscillators is established (thick line in figure (9]^ d)), and backward connection is gradu- 
ally removed (thin line in figure [9t^d)), similar to figures [5t^c, d). As a result, frequency 
synchrony is reached (figure [9]^f)). Based on figures [9](e) and[9]^f), which respectively 
show L and r for different g{0), we estimate Qc-stdp — 7, which is much smaller than 
= 38. Note that g{0) = 7 is a marginal case, which yields a long transient before 
frequency synchrony is reached. Frequency synchrony is not induced with, for exam- 
ple, g{0) = 10 < gc if synapses are frozen. These numerical simulations confirm that 
the results derived in the previous sections apply to networks of pulse-coupled spiking 
neurons. 

6 Discussion 

We have shown that asymmetric STDP greatly reduces the threshold for frequency 
synchronization of neural networks with a pacemaker. This reduction is efficient 
particularly when LTP and LTD are nearly balanced, as assumed for stabilization 



of synaptic weights in previous literature (Kempter et al., 1999 Song et al., 2000 



Song and Abbott, 2001[ |van Rossum et al., 2000]). Our analytical results for two 



oscillator networks provide theoretical understanding of STDP-induced synchrony of 



two-body networks with real neurons (Nowotny et al., 2003) and with Ho dgkin- Huxley 



neurons ( Zhigulin et al. , 2003 ) . Our numerical results for large networks extend ear- 



lier numerical simulations ( Zhigulin and Rabinovich, 2004 ) and provide mechanisms of 
synchrony. 

More microscopically, we have shown that STDP guides formation of feedforward 
networks originating from the pacemaker (figure [5]^g) ) . By eliminating backward con- 
nection, frequency synchrony is promoted in terms of required synaptic weights. Net- 
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works self-organize by asymmetric STDP so that upstream neurons entrain downstream 
neurons. Even though engineered learning algorithms can promote formation of feed- 



forward networks ( Kori and Mikhailov, 2006 ), asymmetric STDP naturally achieves 
this goal. Facilitation of frequency synchrony does not occur for symmetric STDP, 
which cannot suppress backward synapses without sacrificing forward synapses. 

In recurrent networks, synaptic delay may destabilize otherwise stable synchrony, 
leading to oscillatory or chaotic population dynamics (e.g. Gerstner and van Hemmen, 
1993; Gerstner, 2000; Timme et al., 2002). However, our numerical simulations suggest 
that this is not the case in our system, which can be explained as follows. In the phase 
oscillator model, the effect of delay can be replaced in a good approximation by the 



phase shift in the coupling function ( |Kori and Kuramoto, 200H|Kori, 2003D . Dynamics 
of the oscillator system under consideration do not change qualitatively for a large class 



of the coupling function ( |Kori and Mikhailov, 2006| ), which effectively includes the case 
of synaptic delay. Although a synaptic delay enlarges the phase difference between 
connected neurons, the oscillator dynamics in our model are thus robust against delay. 
Another possible complication is that synaptic delay may change synaptic evolution 
because it could interact with the learning window. However, since the delay simply 
increases the phase difference between connected neurons, the causality of spike timing 
does not change even with delay, as corroborated by our numerical experiments. 

In terms of network structure, the feedforward structure is distinct from pruning 
of synapses in a predefined unidirectional network with many presynaptic neurons 



projecting to a single postsynaptic neuron (Kempter et al., 1999: Song et al., 2000 



Song and Abbott, 2001 |van Rossum et al., 2000[). It also differs from multipar- 



tite networks each part of which forms a cluster of synchronously firing neurons 



(Horn et al., 2000 Levy et al., 2001 ). In a sense, feedforward structure and hierar- 



chy are straightforward consequences of asymmetric STDP ( Song and Abbott, 2001 ), 
which opts for causality. We stress that feedforward structure is naturally organized 
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when a network has a distinct pacemaker. The idea of growing feedforward structure 
by unsupervised learning dates back to pioneering work by Bienenstock (1991, 1995), 
which employed Hebbian plasticity. We have analyzed the network formation in detail 
under asymmetric STDP. In real neural networks, there may be multiple pacemakers 
as well as a huge number of follower neurons. It is straightforward to extend our results 
to the case in which a collection of neurons in a network serves as a pacemaker. Pace- 
maker neurons are relevant to, for example, regulation of respiration, internal clock, 
and Parkisonian diseases (see section [1] for references). In these brain regions, pace- 
maker neurons may recruit downstream neurons for frequency synchrony in order to, 
for example, amplify rhythmic activity. 

Formation of feedforward structure could occur even when predetermined pacemak- 
ers are absent. In this case, neurons with relatively high natural frequencies may play 
a role of pacemaker. Backward connection to these fast neurons, which would per- 
turb their periodic firing, can be eventually eliminated by asymmetric STDP. Then, 
the fast neurons can serve as distinct pacemakers. Regardless of the initial presence 
of pacemakers, asymmetric STDP creates frequency synchrony, which can be called 
feedforward synchrony. This mechanism of synchrony differs from that of synchrony 
based on mutual coupling ( [Kuramoto, 1984p . 

Our results do not suggest that asymmetric STDP promotes phase synchrony, 
namely, spike synchrony. This is in contrast to the finding that phase synchrony is 



caused by asymmetric STDP ( Karbowski and Ermentrout, 2002 ). In their work, fixed 
inhibitory coupling as well as excitatory coupling with asymmetric STDP was assumed. 
Perfectly balanced LTP and LTD, at least as the average, is a key condition for the 
maintenance of bidirectional connectivity and phase synchrony. By contrast, we as- 
sumed that LTD is stronger than LTP. This yields a considerable decrease in the 
threshold for synchrony. However, this synchrony is frequency synchrony but not spike 
synchrony. 
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With symmetric STDP, neurons whose spike times are close are hkely to bind 
together. Then, in addition to frequency synchrony, approximate phase synchrony 
whose time resolution is specified by the width of the learning window can develop 



(Seliger et al., 2002). In some situations, neurons divide into clusters in each of which 



rough spike synchrony is maintained ( Masuda and Aihara, 2004 ). Unlike asymmetric 
STDP, symmetric STDP does not lessen the threshold for synchrony. 

When feedforward frequency synchrony is achieved, neurons at different distances 
from the pacemaker fire asynchronously. On top of that, phase synchrony can be ob- 
served for neurons receiving the common signal. For example, the neurons directly 
connected to the pacemaker are excited by the common drive from the pacemaker, 
so that they are synchronized in phase. Likewise, neurons with the same distance 
from the pacemaker tend to fire simultaneously. Indeed, figure [6](b) and its mag- 
nification in figure |6]^c) indicate that clusters of phase-synchronized neurons can be 



aligned according to the distance from the pacemaker (Kori and Mikhailov, 2004). 
Even though spike time difference between neurons with different distances is usu- 
ally small, the order of firing is fixed and reproducible. The pacemaker triggers a 
volley of spikes, which travels down the hierarchy delineated by the distance. This 
phenomenon is consistent with propagation of synfire volley through feedforward 
neural networks in the excitable regime (Bienenstock, 1995 Diesmann et al., 1999 



Reyes, 2003t [Vogels and Abbott, 2005D. However, stably embedding synfire volley in 



recurrent networks is usually difficult ( [Mehring et al., 2003D . It needs, for example, se- 
lective enhancement of forward synapses by 10-fold, which corresponds to large evoked 
EPSPs of 8 mV ( Vogels and Abbott, 2005 ). With asymmetric STDP, forward synapses 
are enhanced. In addition, automatic elimination of backward synapses appreciably 
lessens the forward synaptic strength (or the size of EPSP) needed for stable synfire 
volley. 

We have shown that the facilitation of frequency synchrony by STDP is robust 
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against some heterogeneity in the inherent firing frequency of the neurons. If synap- 
tic delays or neurons are strongly heterogeneous, we would obtain more complex but 
reproducible spatiotemporal spike patterns (Izhikevich et al., 2004 Izhikevich, 2006 



Lengyel et al., 2005). 



Oscillatory neurons can model, for example, temporal coding of place cells 



( Mehta et al., 2002 ) and hippocampal associative memory ( Lengyel et al., 2005 ). By 
contrast, many neural circuits operate in the excitable regime, in which neurons 
are not spontaneously oscillatory. Investigation of the excitable case is warranted 
for future studies. However, we believe that the conclusion that asymmetric STDP 
but not symmetric STDP induces feedforward synchrony generalizes to the ex- 



citable case, as is consistent with previous numerical work ( Song and Abbott, 2001 



Zhigulin and Rabinovich, 2004). 
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Figure 1: Asymmetric (solid line) and symmetric (dashed line) learning windows of 
STDP as a function of tpost — tpre, namely, the postsynaptic spike time relative to the 
presynaptic spike time. 



(a) 




Q CO 
pacemaker oscillator 




oscillator oscillator 



Figure 2: Schematic diagrams showing (a) the network of one pacemaker and one 
oscillator, and (b) the network of two oscillators. 
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Figure 3: Qcstdp for the network with one pacemaker and one oscillator. The evaluation 
by equation f lTB]) (solid line), that by equation f lTS]) (dotted line), and Qcstdp determined 
by numerical simulations of the coupled phase oscillators (circles) are compared. We 

set Qmax — 1-25. 




Figure 4: The critical 5^2 (0) as a function of gi{0) for the network with two oscillators 
(and no pacemaker). Three lines correspond to A'^/A" — 0.96 (thickest hne), 0.9, and 
0.8 (thinnest line). 
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Figure 5: Results for 100 randomly coupled oscillators subject to asymmetric STDP. 
Evolution of the synaptic- weight order parameters are shown for (a, b) ^'(0) = 0.7 
and (c, d) g{0) = 1.5. In (a) and (c), Gf (thick sohd lines), Gb (thin solid lines), Gi 
(moderate solid lines), and the average weight (dotted lines) are shown. In (b) and (d), 
(thick lines) and G^ (thin lines) are indicated. Time courses of (e) L and (f) r are 
compared for g{0) = 0.7, 0.9, 1, 1.2, and 1.5. In (e), lower lines correspond to larger 
g{0). (g) Final network structure for g{0) = 1.5. Only the synapses e E with 
gji > g{0) are presented. The pacemaker is labeled P. Forward edges and backward 
edges are indicated by thin lines and thick lines, respectively. The network is drawn 
by Pajek (fhttp: / / vlado.fmf .uni-lj .si/pub /networks/pajek / ) . (h) Time courses of r for 
some values of g{0) when the oscillators are heterogeneous. 
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Figure 6: Rastergrams of the oscillators under asymmetric STDP in (a) initial and 
(b) final cycles. We set g{0) = 1.5. The oscillators are aligned according to their 
distances li from the pacemaker, which is calculated at time in (a) and 19995 in (b). 
After sufficient time, li is quantized, and the values of li are shown in (b). (c) is a 
magnification of (b). 




Figure 7: Results for 100 randomly coupled oscillators subject to symmetric STDP. 
Evolution of the synaptic weights are shown for g{0) = 140 (a, b) and g{0) = 150 (c, 
d). See the caption of figure O for legends. In (c), Gf (thick sohd line) and Gb (thin 
solid line) overlap almost completely. Time courses of (e) L and (f) r are compared for 
g{0) = 140, 145, 146, 148, and 150. In (e), lower lines correspond to larger 5^(0). (g) 
Final network structure for g{0) = 150. 



-30 

> 




49.5 49.6 49.7 49.8 49.9 50 
time (sec) 



Figure 8: Sample traces of Vq (solid line) and Vi (dashed line) when the spiking neurons 
are uncoupled. 
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Figure 9: Results for 100 randomly coupled spiking neurons subject to asymmetric 
STDP. Evolution of the synaptic weights are shown for (a, b) g{0) = 5 and (c, d) 
g{0) = 10. See the caption of figure [5] for legends. Time courses of (e) L and (f) r are 
compared for g{0) = 5, 6, 7, 8, and 10. In (e), lower lines correspond to larger g{0). 



